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ABSTRACT 

We present an analysis of the starspots on the active M4 dwarf GJ 1243, using four years of time series 
photometry from Kepler. A rapid P = 0.592596 ± 0.00021 day rotation period is measured due to 
the ~2.2% starspot-induced flux modulations in the light curve. We first use a light curve modeling 
approach, using a Monte Garlo Markov Chain sampler to solve for the longitudes and radii of the 
two spots within 5-day windows of data. Within each window of time the starspots are assumed to 
be unchanging. Only a weak constraint on the starspot latitudes can be implied from our modeling. 
The primary spot is found to be very stable over many years. A secondary spot feature is present 
in three portions of the light curve, decays on 100-500 day timescales, and moves in longitude over 
time. We interpret this longitude shearing as the signature of differential rotation. Using our models 
we measure an average shear between the starspots of 0.0047 rad day“^, which corresponds to a 
differential rotation rate of AH = 0.012 ± 0.002 rad day“^. We also fit this starspot phase evolution 
using a series of bivariate Gaussian functions, which provides a consistent shear measurement. This is 
among the slowest differential rotation shear measurements yet measured for a star in this temperature 
regime, and provides an important constraint for dynamo models of low mass stars. 


1. INTRODUCTION 

For low-mass, fully convective stars, the nature of the 
magnetic dynamo and the role of differential rotation is 
not so clear. Some radial and surface differential rota¬ 
tion is expected to exist, due to the combination of ro¬ 
tation and convection. However, despite the deep con¬ 
vective zones of M dwarfs, their long convective turnover 
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, (the H etiectj. S 

mce these stars 


are nearly or tul 
therefore lack a “tachocline” interface region in which 
to store toroidal magnetic field, the dynamo mechanism 
must be fundamentally differen t than the popular afl dy¬ 


namo model for the Sun (e.g. Parker] 1955 Schrijver & 
Zwaan|2000 ). Instead, this convectively driven process is 


known as an dynamo. For rapidly rotating M dwarfs, 
the magnetic field strength is expected to be increased, 
which suppresses di fferential rotatio n and forces nearly 


solid-body rotation (Browning 20081. Without strong ra¬ 
dial or surface differential rotation to organize the global 
magnetic field, activity cycles may not be present, and 
the surface magnetic topology is predicted by some mod- 
els to be highly non-axi symmetric and multipolar (e.g. 


Chabrier fc Kiik^|2006 1. 


However, observations of many low mass stars reveal 


ing low-mass stars show indications of polar “star spot 
caps” possibly due to this large scale dipolar fiel d (Do- 


nati fc Collier Gameron|[l997 Morin et al.||2008'a ), while 
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others do not (|Barnes et al. 2004 Morin et al. 2010). 
Though differential rotation is expected to play a lesser 
role for these rapidly rotating low-mass stars, even small 
amounts of differential rotation may help to organize the 
chaotic, driv en magnetic fields in to a cohe rent, ax- 


isymmetric field (Kitchatinov & Olemskoy 20111, capable 
of producing very long-lived polar spot features. There¬ 
fore, given the wide variety of observed surface magnetic 
topologies, and the complex inter-dependence of rota¬ 
tion, differential rotation, and the magnetic field, mea¬ 
suring differential rotation rates for low mass stars is a 
high priority for constraining dynamo theory. 

Rotation can now be measured with relative ease for 
many stars, for example using spectral line broadening 
that produces v sin i measurements, or periodic flux mod¬ 
ulations due to starspots in precision space-based time se¬ 
ries photom etry. Data from the Kepler mission (Borucki 


et al. 2010) has revolutionized the study of stellar ro- 


tation using starspot modulations, with tens of t hou¬ 
sands of stars having measuring rotation periods (Mc¬ 


Quillan et al. 2014 Reinhold et al. 2013), and has re- 


vealed starspot properties tor stars ranging from solar 


et al 


mass (Bonomo & Lanza 2012) to brown dwarfs (Gizis 
tierentic 


Morin et al. 

2008bD and prominent long-lived starspot 

features (e.g 

Barnes et al. 

2005 

). Some rapidly rotat- 


Ditierential rotation, however, is notoriously difficult 
to detect for stars. Spectral techniques can trace active 
regions at different latit udes for stars with lo wer activity 


levels such as the Sun (Bertello et al. [2012 ). Detecting 
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tween subsequent visits. Photometric surveys may be 
able to produce diff erential rotation rates for an ensem¬ 
ble of active stars (Reinhold & Reiners 2013). A re¬ 
cent blind survey of competing techniques tor detecting 
rotation and differential rotation from model photom¬ 
etry showed excellent agreement in recovering rotation 
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periods froni active stars. However, a complex degener¬ 
acy was found between differential rotation rate, starspot 
lifetimes, and the number of starspots present, and little 
agreement b etween competing me thods and the model 
light curves ( Aigrain fc et al.|20T5 ). In general, methods 
for detecting differential rotation in photometry follow 
one of two approaches: 1) Fourier methods, which mea¬ 
sure the broadening or splitting of peaks in the power 
spectrum, auto-correlation function, or periodogram, or 
equivalently b y decomposition of t he light curve using 
sine functions (Reinhold et al.|2013|). These methods uti¬ 
lize the entire light curve at oiice, and are efficient for an¬ 
alyzing large volumes of data from many stars, but may 
suffer more from the degeneracies mentioned above. 2) 
Tracking specific starspot features e ither via light curve 
inversion (Roettenbacher et al .|2013). or light cur ve mod¬ 
eling for individual starspots (jFrasca et al.||2Ull|). These 
methods are more computationally expensive, but their 
results seem robust for rapidly rotating stars with long- 
lived spots. 

In this paper we venture into a relatively new region of 
starspot evolution parameter space, detecting very grad¬ 
ual differential rotation and spot decay for a rapidly ro¬ 
tating M dwarf. The fast time cadence and continuous 
monitoring provided by Kepler, along with a short stel¬ 
lar rotation period, allow us to trace small changes in 
starspot p hase and amplitude over long periods of time. 
In Sectionwe describe our target, the active M dwarf 
GJ 1243, and the previous investigations of this low- 
mass star with Kepler. Our detailed light curve mod¬ 
eling is presented in Section We trace small changes 
in starspot phase over four years, and interpret this as 
a signature of differential rotation in Section A sim¬ 
pler approach to detect this slow differential rotation by 
modeling the phase evolution with Gaussians is given in 
Section]^ We place the differential rotation signal from 
GJ 1243 in the context of other cool stars, and compare 
the Kepler photometric results with older ground-based 
data in Section]^ Finally, in Sectionj^we provide a sum¬ 
mary of our results, and discuss the great potential for 
understanding starspots and the stellar dynamo still to 
be realized from the unique photometric Kepler database 
and future missions. 

2. GJ 1243 

The target of our study is the nearby mid M dwarf, G J 
1243 (Kepler ID ^ 09726699). This star has a short rota¬ 
tion period of 0.5926 days that h as been noted in previ¬ 


ous studies of Kepler light curv es (Savanov Sz Dmitrienko 


2011 

McQuillan et al.| 2013). 'The 

spectral type has 

been measured as M4 (|Hawley et al. 

2014|), placing GJ 

1243 near the tully-convective boundary where stars are 


expected to remain magnetical l y active for many Gyr 


(N. Reid & S. L. Hawley 
ing th e parallax distance o 
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(|2005 ), the apparent AT-band magnitude from Zacharias 
et al. ( 2013 ), and the M/^-mass relation from Deltosse 
et al.|(|2D00l), we estimate a mass for GJ 1243 of 0.24 Mq. 
The convective turnover timescale for stars in this mass 


range (assuming M=O.235M0) from Kiraga & Stepien 
(|2007|) is quite slow at Tconv ~ 70 days. Comparing this 


timescale to the rotation period, we find GJ 1243 has 
a very low Rossby number of Ro = ProtjTconv ~ 0.008. 
Lucky imaging of GJ 1243, as well as ground based spec¬ 


troscopy, have shown no indications of a binary compan¬ 
ion (Wisniewski 2015 in preparation). In addition, GJ 
1243 has been the subject of detailed flare activity stud¬ 
ies with Kepler data, producing the largest ca talog of 
M dwarf flares ever observed for a single star (Hawleyl 


et al.||2014 Davenport et al. 2014). In this section we 
describe our treatment ot the Kepler data for this ac¬ 
tive M dwarf, removing systematic trends from the light 
curve, and detecting a periodic signal. 

2.1. Kepler Long Cadence Data 

The Kepler light curve for GJ 1243 contains dramatic 
stellar variabi l ity in the form of flares and starspots. 


Ramsay et al. (2013) have examined the flare energy dis- 
tribution using one qua rter (Q6) of data f r om K epler. 


Davenport et al. (2014) and Hawley et al. (2014) used 
11 months ot Kepler short cadence (1-minute) data for 
GJ 1243, over 300 days worth in total, to robustly mea¬ 
sure the flare rate and develop a statistical understand¬ 
ing of the flare morphology from this very active dwarf. 
For these flare studies the starspot signature had been 
treated as a noise source to be smoothed out. 

In the present investigation, we utilized all available 
long cadence (30-minute) Kepler data for GJ 1243 to 
study the evolution of the starspots while minimizing the 
impact of small amplitude flares. GJ 1243 was observed 
in 14 separate quarters of Kepler data (Q0-Q6, Q8-Q10, 
Q12-Q14, and Q16-Q17), spanning over four years of 
observation (MJD 54953.04 through 56423.50). We used 
the most recent reduction of the Kepler data available, 
including the “PDG-MAP” Bayesian de-trending analy¬ 
sis from (Smith et al. 2012). The entire 4-year catalog 
PDG-MAP light curve tor GJ 1243 is shown in Figure 
[T] Data from Q7, Qll, and Q15 was not available due 
to the failure of GGD Module 3 in 2010, which GJ 1243 
resided on for one quarter of the year. 

In Figure large discontinuities in the flux are appar¬ 
ent between quarters, as well as systematic trends in the 
mean flux within quarters. These long timescale varia¬ 
tions are systemic to Kepler data, due to spacecraft drift 
and calibration limitations, and are not astrophysical. 
For every quarter, we fit and subtracted low order (lin¬ 
ear or quadratic) polynomials from the data to remove 
these systematic errors and discontinuities. Because the 
stellar rotation period is so short, and each quarter con¬ 
tains on average ^150 rotations, these polynomial fits 
do not affect the starspot signal on the timescales we are 
interested in. 

Large amplitude flares were also present in our data, 
visible as positive flux excursions throughout the light 
curve in Figure While the short cadence Kepler data 
for GJ 1243 is a tr easure trove for flare studies (e.g. Dav- 
enport et al.|2014|) , only the lar gest energy flares are vis ¬ 
ible m the 3U-minute data (see Walkowicz et al. |2011 ). 
To remove the flares from our aiialysis, we smoothed the 
light curve with a 12-hour “boxcar” filter, and then dis¬ 
carded epochs with fluxes that deviated by more than 
0.3% from the smooth flux. This boxcar smoothing was 
only used to remove outlying epochs, and was not used 
in our starspot analysis. These smoothing values were 
arrived at by eye to remove the most dramatic flares 
and outliers in the data. As this was not a compre¬ 
hensive outlier removal scheme, some small amplitude 
flares and data systematics remained in the light curve. 
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Figure 1. The long cadence PDC MAP light curve for GJ 1243. Pixel shade (light to dark) indicates the density of epochs. Breaks in 
the light curve due to quarterly spacecraft rolls are indicated (grey dashed lines). 


These small amplitude excursions occurred stochastically 
throughout the light curve, had no dependence on rota¬ 
tional phase, and therefore did not affect our spot model¬ 
ing results. As discussed in Lurie et al. (2015), saturation 
can affect the Kepler light curves for ilare stars during 
the brightest flare events. However, the starspot modu¬ 
lations for GJ 1243 were low amplitude, and the quies¬ 
cent flux level was not near the saturation limit. While 
the brightest flux excursions due to flares may be af¬ 
fected by saturation, our starspot analysis is not. Our fi¬ 
nal, inter- and intra-quarter polynomial detrended, flare- 
cleaned light curve for GJ 1243 contained 47,478 epochs 
of data over the four years of Kepler long cadence obser¬ 
vations. 


2.2. Periodic Signal 

The rapid rotation of GJ 1243 was first de tected from 
period ic flux modulations due to starspots by [Irwin et al.| 
(2011) u sing ground-based photometry from the MEarth 
proie ct (Nutzman & Charbonneau 2008 Irwin et al. 
2009). lolfowing the initial QU release oi Kepler data, 
Savanov & Dmitrienko (2011) published the first analysis 
ot the starspots on GJ 1243, using 44 days of continu¬ 
ous long cadence data. They reported a rotation period 
of 0.593 days for GJ 1243, and found that GJ 1243 ex¬ 
hibited two starspot features, separated in longitude by 
203°, and both stable in position over the 44 days of ob¬ 
servation (equal to ~74 rotation periods). The starspots 
covered 3.2% of the visible stellar surface, with a modest 
amount of growth reported over QO. 

A study searching for rotation periods using the auto¬ 
correlation fun ction for ^2500 Kepler M dwarf stars was 
carried out by McQuillan et al. (2013) using 10 months 
of Kepler photometry. They reported a rotation period 
of 0.593 days for GJ 1243 as well. However, a larger 
scale analy s is of o ver 40,000 active Kepler stars by [Rein-] 
hold et al. (2013 ), using the Lomb -Scargle periodogram 
method (|Lomb||l976l |Scargle||1982|), did not report a ro¬ 
tation period tor GJ 1243, as the star’s rapid rotation 
was below their period cutoff. 

These previous studies of GJ 1243 only reported the 
stellar rotation period to an accuracy of 0.001 day (^86 
seconds). With such a short rotation period for GJ 1243, 
an error of 0.001 days would result in phase-folded data 
being out of phase by an entire rotation within one year. 


Thus, to measure any real phase evolution of the starspot 
features over four years we must determine the most ac¬ 
curate mean rotation period possible. We computed the 
normalized Lomb-Scargle periodogram using the entire 
4-year detrended long cadence light curve, using no fre¬ 
quency oversampling or smoothing. The strongest peak 
in the resulting periodogram was very narrow, and had a 
period of 0.592596 days. We then computed the Lomb- 
Scargle periodogram over each of the 14 quarters of data 
individually. The mean period from all quarters we re¬ 
covered was 0.592673 days, which was only ^6.5 seconds 
longer than the period found from all quarters simul¬ 
taneously. These 14 period estimates had a standard 
deviation of 0.00021 days, or about 18 seconds, which 
we adopt as the period uncertainty. Since the rotation 
period is very stable over the course of the Kepler ob¬ 
servations, we assume the period determined from the 
entire light curve, P = 0.592596 ± 0.00021 days, for our 
analysis. 

We then empirically defined the ephemeris of the 
flux minimum by phase-folding the entire Kepler light 
curve at this rotation period. The phase of flux min¬ 
imum was fit using a least squares regression with a 
Gaussian function, which determined an ephemeris of 
to = 2454833.11567807 ± 0.00015. In Figure we 
show median-smoothed 10-day windows of the entire GJ 
1243 light curve, phased using this rotation period and 
ephemeris. The primary dip in brightness stays fixed 
near Phase=0 over the 4 years of observation, which is 
due to the primary starspot. Slow evolution in both 
phase and amplitude of the secondary starspot feature 
is clearly seen. The secondary starspot is almost entirely 
absent at Time ~ 700 days (using units of time as BJD 

- 2454833.11567 days), while the primary and secondary 
starspots appear to have nearly equal amplitudes at Time 

- 1100 . 

Given the long starspot evolution timescale and short 
rotation period for GJ 1243, along with the nearly contin¬ 
uous Kepler light curve for most of the 4-year timespan, 
we are able to study the change in starspot properties in 
much higher temporal detail than illustrated in Figure 
Using 10-day windows of time, we show the 4 year con¬ 
tinuous phase evolution of flux from GJ 1243 in Figure 
1^ For visual clarity the data is folded twice in phase. 
White vertical gaps correspond to quarters with no Ke- 
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Figure 2. Phase-folded, median smoothed light curves for GJ 
1243 from 10-day windows of time, showing the slow evolution of 
the starspot modulations over time. The vertical position for each 
curve corresponds to the start time of the 10-day window on the 
left axis. Each time window is scaled to the same relative flux, 
shown on the right axis. The primary dip, centered at Phase=0 
corresponds to the long-lived starspot. 


pier data, as seen in Figure Each column of pixels in 
this phase versus time flux map contains data spanning 
10 days. This binning resulted in ~16 rotation periods 
per column, with an average of over 400 data points. 
Each row spans 0.04 in Phase, or equivalently a 14.4° 
slice in longitude. The median value for the flux within 
each (time, phase) pixel corresponds to the shading, with 
the darkest regions corresponding to a flux 1.5% below 
the median value, and the lightest pixels 1.5% above the 
median flux. 

The dark band centered at Phase=0 in Figure]^ which 
extends throughout the timespan of the data, is due to 
the primary starspot. This feature does not significantly 
change in phase over the course of our data. The flux 
amplitude for the primary spot is also nearly constant. 
There is an apparent change in the starspot flux ampli¬ 


tude around day ~500 and day ~900, due to the the 
presence of the secondary starspot combined with the 
systematic errors in the flux calibration. The starspot 
features seen in Figure|^are very large compared to spots 
seen on the Sun, appearing to span 50-90° in longitude. 
The detailed geometries of these features cannot be de¬ 
termined from this phase versus time flux map, and each 
observed “starspot” may in fact be a large spot group. 
Additionally, we cannot constrain the total starspot cov¬ 
erage, which may include many smaller spots and active 
regions across the entire stellar surface. Instead we are 
observing the total flux asymmetry due to these spots or 
spot groups. 

The secondary starspot feature continuously changes in 
both phase and flux amplitude (equivalently pixel shade) 
in this diagram. This secondary feature seems to emerge 
and decay at least three times over the 4-year dataset, 
each time appearing nearly on the opposite hemisphere 
of the star and evolving towards the primary starspot. 
Note, a decrease in phase corresponds to a starspot ad¬ 
vancing in longitude in the direction of rotation over 
time. We interpret the slow, linear phase evolution of 
the secondary starspot to be the signature of differential 
rotation. 


3. MODELING THE LIGHT GURVE 

To quantitatively trace the differential rotation on GJ 
1243, we must determine the precise sizes and positions of 
the starspots over time. To accomplish this we performed 
a detailed fit to the Kepler light curve using the starspot 
modeling software from L. Hebb (2015 in preparation). 
Here we give a brief overview of this light curve modeling 
program, as well as our specific use with the GJ 1243 
system. 

The starspot modeling code simulates the star as a 
sphere with uniform surface brightness and limb darken¬ 
ing onto which circular, gray starspots are fixed. Limb 
darkening is implemented by treating the star as a se¬ 
ries of overlapping, concentric circles with brightness val- 
ues defined by the ^coefficient limb darkening model of 
Glaret & Bloemen (20111. Starspots are modeled as non- 
moving circular regions with a fixed flux contrast relative 
to the photosphere, and may be placed anywhere on the 
stellar surface. At each time step in the input light curve, 
as the model star rotates, the code calculates the flux 
blocked by the spots rotating in and out of view, and 
thus generates a synthetic light curve. 

The program can generate a synthetic light curve for a 
single star, with or without a transiting exoplanet, and 
with the spin-axis of the star and orbital axis of the 
planet in any orientation (aligned or misaligned). To 
derive the properties (latitude, longitude and radius) for 
a number of spots that best reproduces the observed flux 
modulations, a comparison is made between the ob¬ 
served data and a synthetic light curve. The model light 
curve generating engine is wrapped with several types of 
Markov Ghain Monte Garlo (MGMG) samplers, includ¬ 
ing an affine invariant MGMG based on |Foreman-Mackey| 
et al. (2012), which explores the parameter space to find 
the lowest and thus the optimum spot properties. 

The program requires that we choose the number of 
spots on the star a priori, and that the spot distribu¬ 
tion remains static. We only analyze a subset of the Ke¬ 
pler data at any one time, using a “window” to model a 
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Figure 3. Continuous phased light curve map for the entire Kepler long cadence dataset. Pixel shade, from dark to light, indicates the 
median flux in each (time, phase) bin. Vertical white gaps correspond to times with no Kepler data, as in Figure]^ Pixels span 10 days 
in time and 0.04 in phase. The starspots are seen as dark regions in this diagram, which evolve in time from left to right. 


timescale over which we do not expect the spots to evolve. 
By sliding this window over the full length of the light 
curve and running the code many times, we fit the entire 
light curve and determine the evolution of the spots. We 
emphasize that the MCMC runs are done independently, 
generating a unique best-fit spot solution within each 
window. This approach of multiple discrete models over 
time avoids parameterizing the starspot evolution with 


analytic fu nctions as has been done previously (e.g. Kip¬ 


ping 20121, which in turn allows us to track non-lmear 
behavior in the size and position evolution of the spots. 
We refer the reader to L. Hebb (2015 in preparation) for 
a description of the full details and capabilities of this 
program, and briefly describe our specihc use below. 

We split the GJ 1243 light curve into windows with 
5 day durations, or approximately 8.4 rotation periods 
at the Kepler 30-minute cadence. The short rotation 
period, combined with the slow evolution of spot fea¬ 
tures seen in Figure resulted in many stellar rotations 
for each window, minimizing the effect of spurious light 
curve features such as flares or small data gaps. Each 
time window was required to contain at least 100 data 
points, or equivalently ~3.5 rotation periods. Each sub¬ 
sequent time window was advanced by 2.5 days, provid¬ 
ing two independent MCMC solutions for each datum. A 
total of 447 such time windows were used spanning the 
14 quarters of data. 

We assumed a fixed flux contrast value of 0.7 for the 
starspots, which is consistent with contrast values seen 
for spots on active giants, as well as the avera ge con¬ 
trast of the solar umbra (e.g. Berdyugina 2005). Note 
that while resulting spot sizes are directly depeiident on 
the contrast value used in our model, the longitude and 
therefore the differential rotation is not affected. The de¬ 


fault value of 100 annuli was used to compute the limb 
darkening. Based on v sin i measurements from echelle 
spectroscopy of GJ 1243 (Wisniweski 2015 in prepara¬ 
tion) and our measured rotation period, we used a fixed 
inclination of 32 degrees. 

For each of the 447 windows of time, we modeled 
the GJ 1243 light curve using two starspots where each 
starspot is defined by fitting three parameters: its lati¬ 
tude, longitude and radius. This was the simplest model 
that was able to reproduce the observed flux modula¬ 
tions for all time windows to high accuracy. We note 
that some time windows were well fit using a single spot 
solution, particularly at Time ^ 800 in Figures and 

where the flux modulation was dominated by a single 
sine-curve like feature. Models with higher numbers of 
spots (three or more) were tested and could easily repro¬ 
duce the observed flux modulations, but were not pre¬ 
ferred when properly compared to the two-spot models 
with fewer free parameters. 

Constraining the latitudes for large starspots is often 
difficult when deriving 2-dimensional starspot conhgu- 
rations from 1-dimensional light curves. There exists a 
well known degeneracy between spot latitude and radius, 
resulting in families of solutions for spots at a given lon¬ 
gitude but a range of latitudes and spot sizes that provide 
equally good fits to the observed light curve. Therefore, 
we chose to fix the latitudes for the two starspots to 
break this degeneracy in our model runs. This does not 
affect our final conclusions because the differential rota¬ 
tion measurements depend only on the derived longitudes 
of the spots. To select latitudes at which to fix the two 
spots in our model, we ran our entire light curve model¬ 
ing analysis for 1/lOth of the time windows, and using 
five configurations of starspot latitudes. For each model 
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configuration one spot was fixed at the steiiar equator 
(0°), and one at a higher iatitude towards the inciined 
pole. The hve higher latitudes spot positions tested were 
(72.8°, 55.6°, 38.4°, 26.9°, 9.8°). Note, given the inclina¬ 
tion of 32°, spots above ~58° would be partially or fully 
visible during the entire stellar rotation, and therefore 
would produce less flux modulation. The starspot longi¬ 
tudes and sizes were allowed to vary in configuration. 

Each resulting set of MCMC solutions produced com¬ 
parably good fits to the light curve, and had the same 
number of free parameters. The individual sine-like mod¬ 
ulations seen in the light curve were not required to cor¬ 
respond to the higher or lower latitude spot in any given 
model configuration. As a result, some models would ex¬ 
hibit a “flip” between spot latitudes for a given feature 
at nearly the same longitude between subsequent time 
windows. This flipping was observed for the two con¬ 
figurations with higher latitude spots (72.8° and 55.6°). 
For our final analysis we chose the solution set with the 
highest latitude configuration that did not exhibit this 
flipping in spot latitudes between subsequent time win¬ 
dows. The two starspots in our analysis were therefore 
fixed at latitudes of 38.4° and 0°. We note our resulting 
longitudinal shear results were insensitive to the latitudes 

chosen. _ 

The affine invariant sampler based on Foreman-Macke^ 


et al. (2012) was employed for each time window, with 
random starting values for the spot radius and longitude, 
but fixed latitudes as described above. For each window 
of data, the MCMC was run for 300 steps using 100 walk¬ 
ers, and the “a scale” parameter was set to 2.0. To carry 
out these independent MCMC realizations efficiently in 
paral lel, we used CONDOR (Litzkow et al.[19M Thain et al. 


20051 to distribute the 447 MCMC explorations across 
180 Linux workstation computer cores. Each window’s 
MCMC chain was converged after 300 steps, and the 
starspot conhguration that produced the best-fit (low¬ 
est x^) solution for each time window was adopted. We 
note that the phase-folded data within each window had 
a scatter about 10 times greater than the typical photo¬ 
metric uncertainty given in the Kepler data. This was 
due to errors in the underlying light curve and limitations 
of our detrending algorithm, as well as small amplitude 
evolution of the starspot features within each window. 
The starspot flux modulation signal was more than 20 
times greater than this scatter. Average reduced val¬ 
ues were ^2 per time window, assuming a 10 times in¬ 
crease in the photometric uncertainty. 

In Figure |4 we show phase-folded light curves for two 
representative time windows of data with their best fit 
solutions overlaid, along with orthographic map projec¬ 
tions of the model stellar surface showing the best-fitting 
spot conhgurations. This map projection demonstrates 
the inclination of the star as well as the relatively large 
size of the starspots. The higher latitude spot appears 
at nearly the same phase (longitude) and size, while the 
lower latitude spot shrinks in radius and advances in lon¬ 
gitude (lower phase) in the direction of the stellar rota¬ 
tion. 


4. QUANTIFYING THE STARSPOT EVOLUTION 

Using the best-fit parameters from each stationary, in¬ 
dependent MCMC model, as in Figure]^ we were able to 
trace the sizes and longitudes of two starspots over the 


entire span of our Kepler data. In Figure!^ we show the 
rotational phase (equivalently the longitude facing the 
observer) for both starspots as a function of time. The 
higher latitude starspot indicated in Figure]^ (orange) is 
very stable in Figure]^ in phase (longitude) with a stan¬ 
dard deviation in longitude of only 4.5% (16 degrees) 
over the four wars of data, and traces the dark band 
seen in Figure centered at Phase = 0. The amplitude 
of this higher latitude spot on the light curve changes 
slowly over the data, with a standard deviation of 34% 
in fractional flux in Figure We refer to this feature as 
the “primary starspot”. 

The “secondary starspot” (purple), however, evolves 
signihcantly in phase across the stellar surface over time 
in Figure!^ This feature corresponds to the lower lati¬ 
tude, equatorial starspot in Figure]^ and traces the tran¬ 
sient secondary features seen in Figure Between Time 
^750 and 900, the two best-fit starspot locations were 
very close in phase, and the variance between solutions 
in subsequent time steps increased for both the primary 
and secondary spots. These correspond to time windows 
where a one-starspot model would be preferred. 

We manually identified two regions in Figurej^that dis¬ 
played nearly constant linear evolution in the secondary 
starspot longitude: Time = 510-630, and Time = 945- 
1400. We interpret these to be the signatures of differen¬ 
tial rotation, with secular spot motions in time. Within 
these time windows we used a non-linear least squares 
first order polynomial fit to measure the linear slopes. 
Lines of best fit for these two regions are shown in Fig¬ 
ure as dashed and solid black lines, and had slopes of 
-0.000927 and -0.000569, respectively. The occurrence of 
these secondary starspots at multiple times within our 
data may in fact be due to a single lower latitude feature 
lapping the primary starspot, but we note the slopes and 
separations in these features in Figure are not consis¬ 
tent with a single spot at a fixed rate of differential ro¬ 
tation. 

The measured slopes were in units of phase day“ ^, and 
corresponded to a rotation shear of AD = 2'K/tiap = 
0.0058 and 0.0036 rad d av~^, using the dehnition from 


Kiiker & Rudiger (2008). Note however this does not in- 
clude any consideration of the starspot latitudes. These 
slopes can also be converted to secondary rotation peri¬ 
ods using the equation 


P, = 


Pn 


I - rriiPo ’ 


( 1 ) 


where nii is the slope of each feature, Pq is the rotation 
period used to phase fold the data to make the figure, and 
Pi is the resulting rotation period. Note by this definition 
a negative slope yields a shorter rotation period. 

Differential rota tion is generally parameterized (e.g. 
Henry et al.||1995 ) as: 


P<l> = Peq/il -ksirr^ cj )), 


( 2 ) 


where is the rotation period at a given latitude 
(0), Peq is the rotation period at the equator, and 
k = AD/Dgq governs the rate of differential rotation 
as a function of latitude. Our model results indicate 
that the period used to phase fold the data in Figure 
corresponds to the higher latitude (38.4°) starspot, 
and assumes the secondary starspot features are on 
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t = 508.1 (days) 


t = 815.6 (days) 





Figure 4. Top: orthographic projections of the model star, with an inclination of 32°, and the best-fit positions for two circular spots for 
the 5 day time window starting at BJD - 2454833.11567 = 508.1 days (left) and 815.6 (right). The direction of stellar rotation is indicated 
by the black arrow. Bottom: phase-folded light curve for the data in the same 5 day time windows, with the best-fit two-spot models 
overlaid (blue solid line), and the contributions from both the higher latitude (orange dashed line) and equatorial (purple dashed line) 
starspots offset for clarity. 
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Figure 5. Top: Continuous phased light curve map, as in Figure]^ with the best-fit solutions from our two spot model overlaid. The 
higher latitude spot shown in Figure (orange open circles) remains nearly constant in phase, while the secondary lower latitude spot 
(purple filled circles) evolves significantly. Linear fits to the phase evolution for the secondary spot are overlaid (black solid and dashed 
lines), which we interpret as differential rotation. Bottom: Fractional flux amplitude of each starspot as a function of time for the best-fit 
solutions from our two spot model. Colors are the same as above. 


the stellar equator. Using an average slope from Fig- 
ure of m = —0.000748, and the phase-folding period 
from §2.2, we estimated an equatorial rotation period 
of Peg = 0.5923336 days via Equation 1. We then 
solved for the unitless differential rotation parameter us¬ 
ing Equation 2, finding k = 0.00114, which corresponds 
to AD = 0.012 ± 0.002 rad day“^. The uncertainty we 
quote here is propagated from the errors in the linear 
least squares fits in Figure 

Assuming the primary spot for GJ 1243 is indeed at a 
higher latitude than the faster rotating secondary spots, 


this behavior is consistent with Solar-like surface differ¬ 
ential rotation where the equator rotates faster than the 
poles. If the primary and secondary starspots are well 
separated in latitude as our model indicates, such a low 
value of shear indicates very weak differential rotation, 
with the star rotating nearly as a solid body. For com¬ 
parison, the Sun’s surface differentia l rotation is much 
stronger, with AD = 0.055 rad day“^ (Berdyugina 2005). 


5. FITTING WITH GAUSSIANS 
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Figure 6. Bivariate Gaussian models of starspot evolution in phase and time (open contours), overlaid on the time-phase flux map from 
Figure]^ A total of four Gaussians defined by Eqn|^were fit, representing one primary and three secondary starspot features. For visual 
clarity we have offset the Gaussian that correspond^o 


correspondsto the primary starspot by 1 phase. 


In the previous two sections we have focused on mea¬ 
suring starspot evolution using a series of sophisticated 
stationary models, and finding the differential rotation 
rate by comparing the position of spots in subsequent 
model realizations. In this section we explore an alterna¬ 
tive method of explicitly determining the starspot time 
evolution, and thus the differential rotation rate, using 
Gaussian functions. 

Rather than modeling the entire light curve directly 
to infer the starspot sizes and positions, as in Figure |4j 
we analyzed the three-dimensional “surface” shown in 
Figures and which traces the flux as a function of 
both time and rotational phase. The data were binned 
in both time and phase (longitude), using bin sizes of 10 
days and 14.4 degrees, respectively. To model this flux 
map we used 2-dimensional bivariate Gaussian distribu¬ 
tions of the form: 


uously wrap from 360 back to 0 degrees. This definition 
enables long-lived starspots with large rates of differen¬ 
tially rotation to “lap” the stellar surface multiple times. 
Each starspot’s evolution is defined by evaluating Equa¬ 
tion over the entire time span of our data, and the full 
range in longitude. The entire flux map in Figure is 
reproduced by summing the 2-D Gaussian functions. 

For our Gaussian analysis we used the data spanning 
from the beginning of the available Kepler data (QO) 
through Quarter 14 (Time ^ 1370). This time range 
was chosen to focus our analysis on the secondary spot 
evolution and differential rotation signal measured above. 
We discarded data from Q16 and Q17 data which showed 
no sign of the secondary starspot. 

We solved for the positions and evolution of four 
starspots over the entire duration of the data, using the 
Python MC MC sampler emcee from |Foreman-Macke}^ 


t' = {t — to) cos 9 — {I — Iq) sin0 
I' = {t — to) sin 9 — (I — Iq) cos 9 


et al. (2012| to explore parameter space for all four si- 
ily. 


(3) 


where F{t, 1) is the flux as a function of time t and lon¬ 
gitude I, A is the flux amplitude of the starspot, r is 
the lifetime of the starspot, b is the scale width of the 
large starspot in longitude, to is the center time of the 
starspot, lo is the center longitude of the starspot, and 
9 is the slope of the spot evolution in units of degrees 
day“^. Here the longitude is a circular coordinate, with 
range between 0 and 360 degrees, and defined to contin¬ 


multaneously. In this four Gaussian model, we consider 
the largest spot (also with the smallest slope in 9) as the 
primary, and the three secondary spots as independent 
spot features, or repeat occurrences of the secondary spot 
discussed before. A third occurrence of the secondary 
spot feature was needed to account for the small feature 
seen around Phase~0.4 at Time~100 days in Figure 
which was not chosen in our conservative by-eye selection 
above. We used rough values for the parameters to seed 
emcee with. The initial seeds for the primary spot were 
9 = 0, to = 800, r = 1300 days, lo = 75 degrees, and 
b = 15. The three secondary spots were all seeded with 
9 = -0.2 degree day“^, r = 100 days, lo = 200 degrees. 
The secondary starspot center times seeds were set to to 
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= 200, 530, and 1200. 

We then ran emcee with 50 walkers for 2000 steps. 
The best fit model from this parameter space search for 
the three secondary starspots is shown in Figure The 
best fit slope for the primary spot was 9 = -0.0011 deg 
day“^, and for the secondary spots (in time order) was 
9 = -0.20, -0.14, and -0.15 deg day“^. These secondary 
spot shear rates corresponded to 2TT/tiap = All = 0.0036, 
0.0025, and 0.0027 rad day“^, or an average of 0.0029 rad 
day“^, somewhat lower amplitude than measured from 
the linear features in Figure Note again this method 
does not constrain the latitudes of the starspots, and so 
the measurement of shear is only a lower limit on the 
true differential rotation rate. 

This approach assumes a priori that the starspot evo¬ 
lution in both time and longitude can be represented by 
a Gaussian function, meaning the spots may only evolve 
linearly in longitude over time. We note the resulting 
estimate for the differential rotation shear rates for the 
secondary starspots are very similar to the values de¬ 
termined when fitting many time-stationary MCMC in¬ 
stances. The Gaussian modeling approach used flux data 
that was binned in time and phase, greatly reducing the 
number of data points to be fit. This entire MCMC anal¬ 
ysis took only a few minutes to compute using a standard 
Linux workstation. We thus propose this to be an effi¬ 
cient means of estimating the starspot lifetimes and dif¬ 
ferential rotation rates in certain cases, which we discuss 
further in Section below. 

6. DISCUSSION 

We have produced two significant results in this work. 
The first is the identification of a long-lived starspot, 
which we attribute to a higher latitude starspot (possi¬ 
bly due to a spot cap or group) on the rapidly rotating M 
dwarf, GJ 1243. The second is a robust measurement of 
weak differential rotation for this star due to a spot closer 
to the stellar equator. In this section we provide addi¬ 
tional context and discussion of these results, and their 
implications for the magnetic field’s surface topology. 

6.1. A Long Lived Starspot 

To further illustrate the remarkable stability of the 
higher latitude spot on GJ 1243, we retrieved ground- 


based light curves from SuperWASP which 

aredated the 

Kepler mission by ^2 years (Butters et al. 

2010 

). This 


SuperWASP public archive photometry was phased us¬ 
ing the period and ephemeris we determined from our 
Kepler light curve. The phase of flux minimum matches 
between these two datasets to within 1%, indicating this 
large starspot has been stable in longitude for more than 
6 years. The amplitudes of the flux modulations between 
Kepler and SuperWASP are only slightly different, with 
the median Kepler variation of 2.19% (averaged over all 
4 years of data), and for SuperWASP of 2.86%. The Su¬ 
perWASP data was taken in the V-band, which is more 
narrow than the very wide Kepler filter. The V-band 
also is centered at a shorter wavelength than the Ke¬ 
pler filter, which is weighted more towards the R-band. 
As a result, we would expect to find larger flux contrast 
between cool spots and the stellar photosphere in the V- 
band than compared to the Kepler filter. However, since 
these observations were not concurrent we cannot rule 


out small differences in the starspot’s physical size over 
time. 


6.2. Differential Rotation in Cool Stars 

The average starspot shear for GJ 1243 observed in 
this paper of 0.0047 rad day“^ corresponds to a differ¬ 
ential rotation rate of AD = 0.012 ± 0.002 rad day“^ 
(assuming the spot configurations used in our models), 
and is one of only a few such measurements yet ob¬ 
tained for low mass, rapidly rotating, fully convective 
stars. In Figure we place this measurement in the 
context of othe r existing observations of stellar differ- 
ential rotation (|Barnes et al. 2005 Morin et al. 2008a I 
along with the empirical extrapolation to cool stars trorn 
Reiners ( 2006 Go l lier C ameron (20071, and models from 
Kiiker & Rudiger 12011). One of the few other objects 
with a robust differential rotation measurement in this 


regime is the cool, rapidly rotating star, V374 Peg (Morin 


et al. 2008a I. These authors employ Doppler Imaging, a 


completely different technique to our own, to derive a 
value for the surface differential rotation of V374 Peg 
(AH = 0.0063±0.0004 rad day“^) that is similar to that 
of GJ 1243 we have measured. This Doppler Imaging 
method assumes a Solar-like differential rotation profile 
as in Equation 2, and simultaneo usly fits for t he sta rspot 
positions, sizes, and shear rates. Lurie et al. ( 2015|) have 
also estimated the starspot shear tor both rapidly ro¬ 
tating components of the M5-I-M5 binary GJ 1245AB, 
using Kepler data and the 2D Gaussian modeling ap¬ 
proach detailed in our Section]^ finding shear rates that 
are comparable to GJ 1243 and V374 Peg. 

The various methods for constraining differential ro¬ 
tation each have unique limitations in their sensitivity 
and degeneracies. Photometric phase-tracking methods 
such as ou rs, frequency splitting approaches (Reinhold 
et al. 2013), as well as spectroscopic line broad ening tech¬ 
niques (Tie, th e Fourier Transform Method; Reiners & 
Schmitt 2003) have been considered as only providing 
lower limits on the true differential rotation ra te, due 
to a lack of constraint on the starspot latitudes (Barnes 
et al. 2005). Statistical correctio ns for this latitudina l 
uncertainty have been developed (Hall & Henry 1994), 
and our models of GJ 1243 provide a weak constraint 
on the spot latitudes. Together, these varied observa¬ 
tions indicate that rapidly rotating, mid M-dwarf stars 
exhibit differential rotation that is significantly weaker 
than on the Sun (0.055 rad day“^) by up to an order of 
magnitude. 

These observational results of low AH values support 
recent theoretical work in this area. Mean field theory 
models predict that for stars at a fixed temperature sur¬ 
face differential rotation rates decrease with faster ro- 
tation (shorter periods), down to periods of a few days 
(IKiiker & Rudiger||2005). In addition, the total ampli- 
tude of the differential rotation decreases with decreas¬ 
ing temperature, as models indicate lower mass main 
sequence stars (Teff < 6000 K) should exhibit lower 
amounts of surface differen tial rotation than hotter stars 
at a fixed rotation period ( Kiiker fc Rudiger||2011 ). 

Dipolar dynamos with strong magnetic held and 
quenched differential rotation have been reported in 
glob al dynamo models of r apidly rotating low-mass stars 
(e.g. Gastine et al. 2013). Furthermore, recent global 
dynamo modeling efforts have also been successful in 
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Figure 7. Average starspot shear from the two linear fits to the 
MCMC light curve models in Figurel^ assuming the primary spot 
was at a latitude of 38.4° and secondary spot at the equator (blue 
diamond). The blue bar extends to the minimum possible am¬ 
plitude of differential rotation for GJ 1243, assuming the primary 
and secondary s pots are at the pole and equator, respectively. For 
comparison, thejCollier Cameron| (|2007|l observe d fit for cool stars 
fblack dashed liriej, theoretical prediction from|Kuker &: Riidigerl 
J2011|| (red dot-das h and blue solid lines), individual stars trom 
|hSarnes et al.| (|2005ll (black circles), and the estimate d shear rate 
tor VdY 4 J-"eg determined using doppler imaging from |Morin et al.| 
< |2008a| l (purple filled circle) are shown. 


produ cing polar starspots self consistently (Yadav et al. 
20151 Furthermore, in hydrodynamic simulations ot ro^ 


tatmg, non-magnetic, solar-type stars, Browning (20081 
found that convection in the outer zone ot the star redis- 
tributes angular momentum quickly giving rise to solar¬ 
like differential rotation. However, when strong magnetic 
fields are introduced, the field lines act to reduce the dif¬ 
ferential rotation by linking together individual regions 
of the stellar interior. As the strength of the magnetic 
field is increased, as is typically seen for more rapidly 
rotating stars, the differential rotation is suppressed to 
almost negligible values. 


The observations from Lurie et al. (2015) also sup¬ 


port this connection between rotation, magnetic field 
strength, and suppressed differential rotation. For the 
M5-I-M5 binary system GJ 1245AB, the more rapidly 
rotating component (GJ 1245A, P=0.26 days) has a 
slightly higher total chromospheric Ha emission flux, and 
significantly less phase evolution of its starspot modula¬ 
tions compared with GJ 1245B, (P=0.71 days). 

6.3. Magnetic Field Topology 

Dynamo models and observations of stars like GJ 1243 
are in agreement that for rapidly rotating, low-mass 
stars, the magnetic field strength is increased and surface 
differential rotation is suppressed. The detailed surface 
topology of the magnetic fie ld is less certain . The con¬ 
vective dynamo models from Browning (2008), for exam¬ 
ple, indicate that rapid rotation produces strong axisym- 
metric magnetic fields. Ear lier models of convective en¬ 
velopes (|Dobler et al.||2006|) also show net axisymmetric 
magnetic helds. Other models of fully convective stars, 
however, ha ve produced fully non-ax i symmetric fields at 
the s urface (Kiiker & Riidigerl 2005 Chabrier & Kiiker 
20061 ). 


The smooth rotational flux modulation observed is due 
to a local feature on the stellar surface. The magnetic 
activity from this star as t raced by flares in Kepler data 


has also be en well studied (Hawley et al. 2014 Davenport 
et al.|20l4 |. No correlation between rotational phase (or 
equivalently longitude) and the occurrence rate or energy 
emitted from flares has been found. The star may be uni¬ 
formly covered by many smaller active regions and spots 
that would not create observed modulations in the light 
curve. The flares would be a result of these small scale 
multipolar magnetic field structures. Similarly, the stel¬ 
lar pole may be entirely covered by a polar spot “cap”, re¬ 
sulting from strong poloidal magnetic flux geometry due 
to the rapid stellar rotation. This poloidal component of 
the magnetic field could also be slightly misaligned from 
the stellar rotation axis, as seen in other convective stel¬ 


lar and planetary dynamos (e.g. |Christensen et al.||2009 
Hull et al. 2013), resulting in the observed light curve 


asymmetries. 

As we noted in Section 2, the Rossby number for GJ 
1243 is Ro « 0.008. Acco rding to the Z DI observations 
of M dwarfs aggregated in|Gastine et al.| (|2013|), for stars 
with Ro < 0.1 both dipole and multipole helds are pos¬ 
sible. In this low Rossby number regime, they also find 
surface differential rotation should be stronger when the 
magnetic field is multipolar. However, we find for GJ 
1243 a very low rate of differential rotation, and large, 
long-lived starspot modulations. As a result, in this con¬ 
text we predict a highly organized and stable dipolar 
magnetic field geometry. This is in agreement with the 
ZDI observations of the similar star V374 Peg, wh ich has 


a large scale dipolar field and long-lived starspots (Morin 


et al. 


2008a). 


7. SUMMARY 

In this paper, we have presented a classic approach of 
phase tracking starspots in light curves, made new by 
the exquisite photometric monitoring from Kepler. By 
tracing the phase evolution for two starspot regions on 
GJ 1243 we have found the smallest amplitude of dif¬ 
ferential rotation rate ever robustly measured for a cool 
star. This p hase-tracking technique is similar to |Henry| 
et al. (1995), and is sensitive to a comparable amp litude 
differential rotation signal to Morin et al. (|2008a). The 
large starspots, or starspot groups, on GJ 1243 are very 
long-lived, with the primary high-latitude spot found to 
be constantly aligned in phase for over 6 years. The sec¬ 
ondary starspot features evolve on timescales of hundreds 
of days in both phase and amplitude. 

There remain many challenges in modeling the 
starspots using broadband light curves alone. For exam¬ 
ple, we have almost no constraint on the actual latitudes 
of the spots. Modeli ng starspots on stars with tran¬ 
siting exoplanets (e.g. |Sanchis-Ojeda et al.||2013D may 
help break many of these degeneracies. As most tran¬ 
siting systems in the Kepler data are aroun d G dwarf 


stars, more M dwarf systems like Kepler 186 (Quintana 
et al.|2014 ) are needed to better understand the detailed 


starspot characteristics of stars across the main sequence. 

Both the light curve modeling MCMC and 2D Gaus¬ 
sian phase-tracking techniques used to measure surface 
differential rotation in this work are best suited for track¬ 
ing long-lived spots on rapidly rotating stars, as in the GJ 
1243 system where we are able to average over hundreds 
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of stellar rotation periods during a starspot’s lifetime. 
We believe the light curve modeling approach provides 
the most robust estimates for starspot properties, but 
note the 2D Gaussian approach is orders of magnitude 
faster to execute. This methodology could be applied to 
hundreds of rapidly rotating active stars in the Kepler 
dataset f or which rotation periods are already known 
(e.g. see Reinhold et al. 2013 |McQuillan et al. 2014). 
There are ~20 other stars in the Kepler data with M 
dwarf colors and estimated rot ation periods shorter than 
1 day (McQuillan et al. 2013). A cursory look at these 
light curves reveals many with dramatic flare activity and 
sinusoidal starspot modulations, as found on GJ 1243 
and GJ 1245AB. The phase versus time diagrams (as 
in our Figure]^ for GJ 1243) for these other rapidly ro¬ 
tating stars show a diverse set of morphologies, ranging 
from even more stable spots than on GJ 1243, to stars 
with faster shear rates and shorter spot lifetimes. 

Finally, we have introduced an efficient technique for 
empirically tracking starspot evolution in the phase ver¬ 
sus time flux map, by fitting bivariate Gaussians to 
model the spot motion and evolution. Studying the 
phase-evolution of the starspots with Kepler for single 
field stars appears to be feasible for stars with fast ro¬ 
tation periods and long spot lifetimes. The Gaussian- 
fitting method presented here has already been applied 
to the Kepler dat a for the active M 5-I-M5 binary sys¬ 
tem, GJ 1245 AB (Lurie et al. 2015). We have pointed 
out many other rapidly rotating low-mass stars in the 
Kepler archive that may be studied with this technique, 
and hope this work will be the beginning of a larger ob¬ 
servational understanding of surface differential rotation 
in cool stars. 
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